This report summarises the analysis and results for the publication “Cognitive benefits of learning additional languages in old adulthood? Insights from an intensive longitudinal intervention study”.

The dataset

63 participants completed the training, two of whom were excluded, because they suffered a minor stroke halfway through the 30 weeks and did therefore no longer qualify as neurologically healthy.

The dataset therefore includes 61 trajectories of cognitive performance over a period of 30 weeks, assessed on a weekly basis.

See Anonymized (2021) for a detailed description of the cognitive tasks and socio-affective questionnaires.

The tests include:

  • Alertness task
    • Extracted measure: Reaction time (RT)
    • Accuracy almost 100%, hence not used.
  • Divided attention task
    • Extracted measures: Accuracy (i.e. sensitivity) & RT.
  • Verbal fluency task (“Regensburger Wortflüssigkeitstest”)
    • Number of words produced in 1 minute
  • Complex working memory task (Operation Span)
    • Extracted measure: Accuracy
    • RT not extracted, because letters had to be entered, so it depends largely on computer literacy
  • Simple working memory task (2-Back)
    • Extracted measures: Accuracy (= sensitivity) & RT

Exploratory Analyses

ID-Differences between Groups

In order to avoid between-group effects based on a selection bias, groups were compared in terms of age, number of years in education, multilingualism and number of regular activities/hobbies:

# See idDiffs.R for code.

load(file = paste(inputFolder, "id_dfLong.rda", sep = "/"))

# remove standardized variables
id_dfLong2 = id_dfLong %>%
  filter(Measure %in% c("ageZ", "educ_yZ", "BLPZ", "Act_totalZ") )

id_dfLong2 = droplevels(id_dfLong2)

levels(id_dfLong2$Measure) <- c("Age", "Education", "Multilingualism", "Number of Hobbies")
id_dfLong2$Group = relevel(id_dfLong2$Group, ref = "ACTV")
id_dfLong2$Group = relevel(id_dfLong2$Group, ref = "PASV")

id_dfLong2 %>%
  ggplot(aes(x = Group, y = Value, fill = Group)) +
  geom_boxplot() +
  facet_grid(~ Measure, scales = "free") +
  theme_bw() +
  ylab("Z-Score") +
  xlab("")+
  scale_fill_brewer(palette = "Paired", limits = levels(id_dfLong2$Group)) 

# Chi-square test for gender distribution
gend = data.frame(LANG = c(14,14), ACTV = c(7,10), PASV = c(5,11))
row.names(gend) = c("f", "m")

chi = chisq.test(gend)

See 5_publPlots.R for descriptive statistics table.


Result:

  • Groups were comparable in their ID variables.
  • The three training groups did not differ with respect to the ratio of male to female participants χ2(2) = 1.48, p = 0.48.



Variable Reduction

Let’s see if the number of variables can be reduced based on correlations. We will also include the predictors “Wellbeing” and “Training Motivation”, as they are likely to correlate.

# correlation matrix between cognitive variables
originalDat = read.csv(paste(inputFolder, "selectedData.csv", sep = "/"))

thisDat = originalDat %>%
  select(sNBack, sOpSpan, sRWT, rtNBack, rtDivAtt, rtAlert, Wellbeing, TrainingMotivation) %>%
  mutate(rtNBack = log(rtNBack), # logarithmize all RTs
         rtAlert = log(rtAlert),
         rtDivAtt = log(rtDivAtt)
)

# make correlation matrix
m = cor(thisDat, method = "spearman", use = "pairwise") # spearman because non-normal distribution

# plot
plot_ly(x = names(thisDat), y = names(thisDat), z = m, type = "heatmap")

Result:

The strongest correlations are between:

  • Wellbeing + Training Motivation: r = 0.7
  • N-Back Accuracy + Operation Span Accuracy: r = 0.34

It makes sense from a theoretical point of view to concatenate Wellbeing & Training Motivation and N-Back Accuracy & Operation Span Accuracy. The resulting variables aresocioAffect and sWM (Working Memory Accuracy), respectively.



Distribution of Cognitive and Socio-Affective Variables

# Make new df in long format with cognitive and socio-affective variables
thisDat = df %>%
  select(userCode, Group, rtAlert, sDivAtt, rtDivAtt, sWM, rtWM, sRWT, socioAffect) %>%
  tidyr::pivot_longer(rtAlert:socioAffect,
                      names_to = "Variable",
                      values_to = "Score") %>%
  group_by(Variable) %>%
  mutate(outlier.high = Score > quantile(Score, .75, na.rm = TRUE) + 1.5*IQR(Score, na.rm = TRUE),
         outlier.low = Score < quantile(Score, .25, na.rm = TRUE) - 1.5*IQR(Score, na.rm = TRUE)) %>%
  ungroup()

thisDat$outlier.color <- NA
thisDat$outlier.color[thisDat$outlier.high] <- "red"
thisDat$outlier.color[thisDat$outlier.low] <- "steelblue"

# Make Variable a factor
thisDat$Variable = as.factor(thisDat$Variable)

# Relevel Group for coloring
thisDat$Group = relevel(thisDat$Group, ref = "ACTV")
thisDat$Group = relevel(thisDat$Group, ref = "PASV")

# Rename factor levels
levels(thisDat$Variable) = c("Alertness RT [log]", "Divided Attention RT [log]", "Working Memory RT [log]",
                             "Divided Attention Acc.", "Socio-Affect", "Verbal Fluency", "Working Memory Acc.")

# Draw plot
plotList = list()

for (i in 1:length(unique(thisDat$Variable))) {
  plotList[[i]] = thisDat %>%
    filter(Variable == unique(thisDat$Variable)[i]) %>%
    ggplot(aes(x = Group, y = Score)) +
    geom_violin(aes(group = Group, fill = Group),
                trim = FALSE) +  
    geom_boxplot(aes(group = Group), width = 0.15, outlier.shape = NA) +
    geom_point(aes(color = outlier.color), data = function(x) dplyr::filter_(x, ~ outlier.high | outlier.low), position = position_jitter(w = 0.1, h = 0.05), alpha = 0.5) +
    stat_summary(mapping = aes(group = Group, fill = Group), fun = mean, geom = "point", shape = 23, size = 2.5) +
    scale_fill_brewer(palette = "Paired") + #Alternatives: "Dark2", "RdBu"
    scale_x_discrete() +
    theme_bw() +
    theme(legend.position = "none") +
    ylab("Score") +
    ggtitle(unique(thisDat$Variable)[i]) +
    theme(plot.title = element_text(hjust = 0.5, face = "bold")) + 
    xlab("")
}

ggarrange(plotList[[1]], plotList[[2]], plotList[[3]], plotList[[4]], plotList[[5]], plotList[[6]], plotList[[7]],
          nrow = 3, ncol = 3)


N_vect = thisDat %>%
  group_by(Group, Variable) %>%
  filter(is.na(Score) == F) %>%
  count()

sumDF = thisDat %>%
  group_by(Variable, Group) %>%
  summarise(N = 0,
            Mean = mean(Score, na.rm = T),
            Median = median(Score, na.rm = T),
            "Std. Dev" = sd(Score, na.rm = T),
            Min = min(Score, na.rm = T),
            Max = max(Score, na.rm = T)) 

sumDF = full_join(sumDF, N_vect, by = c("Variable", "Group"))

sumDF = sumDF %>%
  dplyr::select(Variable, Group, N = n, Mean, Median, "Std. Dev", Min, Max)

sumDF %>%
  kable(digits = 2, format = "html", booktabs = T,
      caption = "Descriptive Summary of Repeated Measures") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
Descriptive Summary of Repeated Measures
Variable Group N Mean Median Std. Dev Min Max
Alertness RT [log] PASV 453 -1.15 -1.17 0.14 -1.39 -0.29
Alertness RT [log] ACTV 489 -1.17 -1.16 0.12 -1.46 -0.67
Alertness RT [log] LANG 809 -1.17 -1.19 0.14 -1.42 -0.62
Divided Attention RT [log] PASV 452 0.07 0.11 0.26 -0.55 0.53
Divided Attention RT [log] ACTV 489 -0.03 -0.05 0.25 -0.58 0.44
Divided Attention RT [log] LANG 805 -0.02 -0.03 0.26 -0.62 0.56
Working Memory RT [log] PASV 450 -0.36 -0.37 0.25 -0.98 0.33
Working Memory RT [log] ACTV 489 -0.32 -0.36 0.27 -1.02 0.41
Working Memory RT [log] LANG 802 -0.44 -0.46 0.25 -0.98 0.49
Divided Attention Acc. PASV 452 6.27 6.78 1.61 1.49 8.22
Divided Attention Acc. ACTV 489 6.87 7.01 1.29 1.57 8.22
Divided Attention Acc. LANG 796 6.93 7.03 1.34 0.23 8.22
Socio-Affect PASV 440 78.39 81.02 17.19 16.80 100.00
Socio-Affect ACTV 457 76.05 79.88 20.82 11.59 100.00
Socio-Affect LANG 785 74.26 77.73 18.53 8.72 100.00
Verbal Fluency PASV 451 12.11 12.10 3.44 4.10 22.17
Verbal Fluency ACTV 488 12.60 12.19 3.10 4.19 21.10
Verbal Fluency LANG 806 13.01 13.10 3.05 4.19 22.51
Working Memory Acc. PASV 455 0.72 0.74 0.15 0.16 1.00
Working Memory Acc. ACTV 489 0.77 0.78 0.14 0.20 1.00
Working Memory Acc. LANG 805 0.80 0.80 0.12 0.39 1.00


Result:

  • There is a clear ceiling effect in the case of Divided Attention Accuracy. The score will be removed from further analysis.
  • The distribution of socio-affect is skewed to the left, i.e. participants were highly motivated.
  • Note: Reaction times are already log-transformed.




Initial Conditions

To quantify initial conditions, the very first time point was removed as familiarization trial (in which participants also received feedback from the investigators). Results did not change noticeably when this time point was included. A GAM was computed for each subject and task, and the score at Time = 0 was retrieved from each of these GAMs. The mean over all five initial scores was defined as the individual’s initial conditions.

Ideally, groups should perform similarly on the cognitive tasks at the start of their respective trainings. Whether this is the case, will be tested here.

The following are boxplots for the overall GAMM intercepts and intercepts extracted for each task individually indicating significant differences in intercept between groups.

load(file = paste(inputFolder, "id_df.rda", sep = "/"))

thisDat = id_df %>%
  group_by(Group) %>%
  mutate(outlier.high = initLevel > quantile(initLevel, .75, na.rm = TRUE) + 1.5*IQR(initLevel, na.rm = TRUE),
         outlier.low = initLevel < quantile(initLevel, .25, na.rm = TRUE) - 1.5*IQR(initLevel, na.rm = TRUE)) %>%
  ungroup()

thisDat$outlier.color <- NA
thisDat$outlier.color[thisDat$outlier.high] <- "red"
thisDat$outlier.color[thisDat$outlier.low] <- "steelblue"

# Relevel Group for coloring
thisDat$Group = relevel(thisDat$Group, ref = "ACTV")
thisDat$Group = relevel(thisDat$Group, ref = "PASV")

thisDat %>%
  ggplot(aes(x = Group, y = initLevel)) +
  geom_violin(aes(group = Group, fill = Group),
              trim = F) +
  geom_boxplot(aes(group = Group), width = 0.15, outlier.shape = NA) +
  geom_point(aes(color = outlier.color), data = function(x) dplyr::filter_(x, ~ outlier.high | outlier.low), position = position_jitter(w = 0.1, h = 0.05), alpha = 0.5) +
  stat_summary(mapping = aes(group = Group, fill = Group), fun = mean, geom = "point", shape = 23, size = 2.5) +
  scale_fill_brewer(palette = "Paired") + #Alternatives: "Dark2", "RdBu"
  scale_x_discrete() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme_bw() +
  xlab("") +
  ylab("Baseline") +
  ggtitle("Overall Cognitive Baselines")



Predictors of Initial Conditions

load(file = paste(inputFolder, "id_df.rda", sep = "/"))

corrDF = melt(id_df, id.vars = c("userCode", "initLevel"),
              measure.vars = c("age", "educ_y",  "BLP", "Act_total", "socioAffect_init"),
              variable.name = "Measure",
              value.name = "Value")

levels(corrDF$Measure) = c("Age", "Education", "Multilingualism", "Number of Hobbies", "Initial Motivation")

# Make a dot and a line color column depending on the significance level
# Check significance level and - if significant - add respective color to new column

ageR = data.frame(dotCol = rep(ifelse(summary(lm(initLevel ~ age, data = id_df))$coefficients[2,4]< 0.05,
                                     "khaki3", "lightgrey"), 
                     length(unique(id_df$userCode))),
                 lineCol = rep(ifelse(summary(lm(initLevel ~ age, data = id_df))$coefficients[2,4]< 0.05,
                                      "cyan4", "darkgrey"), 
                               length(unique(id_df$userCode))))
educR = data.frame(dotCol = rep(ifelse(summary(lm(initLevel ~ educ_y, data = id_df))$coefficients[2,4]< 0.05,
                                      "khaki3", "lightgrey"), 
                               length(unique(id_df$userCode))),
                  lineCol = rep(ifelse(summary(lm(initLevel ~ educ_y, data = id_df))$coefficients[2,4]< 0.05,
                                       "cyan4", "darkgrey"), 
                                length(unique(id_df$userCode))))
blpR = data.frame(dotCol = rep(ifelse(summary(lm(initLevel ~ BLP, data = id_df))$coefficients[2,4]< 0.05,
                                      "khaki3", "lightgrey"), 
                               length(unique(id_df$userCode))),
                  lineCol = rep(ifelse(summary(lm(initLevel ~ BLP, data = id_df))$coefficients[2,4]< 0.05,
                                       "cyan4", "darkgrey"), 
                                length(unique(id_df$userCode))))
actR = data.frame(dotCol = rep(ifelse(summary(lm(initLevel ~ Act_total, data = id_df))$coefficients[2,4]< 0.05,
                                      "khaki3", "lightgrey"), 
                               length(unique(id_df$userCode))),
                  lineCol = rep(ifelse(summary(lm(initLevel ~ Act_total, data = id_df))$coefficients[2,4]< 0.05,
                                       "cyan4", "darkgrey"), 
                                length(unique(id_df$userCode))))
motivR = data.frame(dotCol = rep(ifelse(summary(lm(initLevel ~ socioAffect_init, data = id_df))$coefficients[2,4]< 0.05,
                                      "khaki3", "lightgrey"), 
                               length(unique(id_df$userCode))),
                  lineCol = rep(ifelse(summary(lm(initLevel ~ socioAffect_init, data = id_df))$coefficients[2,4]< 0.05,
                                       "cyan4", "darkgrey"), 
                                length(unique(id_df$userCode))))

tmp = rbind.data.frame(ageR, educR, blpR, actR, motivR)
corrDF = cbind(corrDF, tmp)

# clean up
remove(ageR, educR, motivR, actR, blpR, tmp)

# Plot

corrDF %>%
  ggplot(aes(x = Value, y = initLevel)) +
  geom_point(aes(color = dotCol)) +
  geom_smooth(method = "lm", alpha = 0.2, aes(color = lineCol)) +
  facet_grid(~ Measure, scales = "free") +
  theme_bw() + 
  ylab("Cognitive Baseline") +
  xlab("") +
  #stat_cor(method = "pearson", 
  #         label.x = 0.25, label.y = 0.9, label.sep = "\n")  +
  scale_color_manual(values = c("cyan4", "darkgrey", "khaki3", "lightgrey")) +
  guides(color = F)

Even though there was no significant difference in ID variables between groups, education (r = 0.35, p = 0.01) and multilingualism (r = 0.41, p = 0) predicted the initial cognitive performance and likely caused the group differences.

The correlation between multilingualism and education was negligible (r = 0.17, p = 0.19)



Differences in Socio-Affect

To ensure that socio-affect was comparable between groups, a GAMM was computed with socio-affect as the dependent variable and individual smooths for each group (as ordered factors). Due to the approaching ceiling effect, the model’s residuals were not distributed normally even after using the scaled-t family:


load(paste(outputFolder, "tables_and_summaries", "socioAffect_OFgamm.rda", sep = "/")) # Model is run in idDiffs.R
# model is called "m1a"

qqp(resid(m1a), main = "Model Residuals Socio-Affect GAMM", ylab = "Residuals", id = F)


The mean socio-affect per group and the respective GAMM smooths are visualized in the following plot. See 2_idDiffs.R for plot code.



Result: All groups were highly motivated, which is why there was no normal distribution and the model could not fit the data perfectly. The model (albeit with non-normal residuals) could not find any significant difference between groups in terms of their initial socio-affect or the development thereof:


# For respective model, run this section:
load(paste(outputFolder, "tables_and_summaries", "summ_socioAffect_OFgamm.rda", sep = "/")) # Model is run in idDiffs.R
summ_m1a
## 
## Family: Scaled t(3,3.859) 
## Link function: identity 
## 
## Formula:
## socioAffect ~ isACTV0 + isPASV0 + s(Time) + s(Time, by = isACTV0) + 
##     s(Time, by = isPASV0) + s(Time, userCode, bs = "fs", m = 1)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   79.181      2.684  29.499   <2e-16 ***
## isACTV01       1.957      4.370   0.448    0.654    
## isPASV01       3.355      4.711   0.712    0.476    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf  Ref.df       F p-value    
## s(Time)            1.000   1.001   1.128   0.288    
## s(Time):isACTV01   1.005   1.006   1.320   0.249    
## s(Time):isPASV01   1.578   1.699   0.587   0.424    
## s(Time,userCode) 526.322 546.000 213.933  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.695   Deviance explained = 68.2%
## fREML =  17913  Scale est. = 1         n = 8490
# Get mean socio-affect per subject
thisDat = df %>%
  group_by(Group, userCode) %>%
  summarise(socioAffect = mean(socioAffect, na.rm = T))

test = kruskal.test(socioAffect ~ Group, data = thisDat)

A Kruskal-Wallis Test did not show any difference in mean motivation between groups χ2(2) = 0.53, p = 0.77.



Summary of Exploratory Analysis

  • Groups were comparable in terms of all ID variables.
  • Group LANG slightly outperformed group PASV at baseline.
  • Baseline differences were predicted by education and multilingualism.
  • All groups were highly motivated and reported high wellbeing.



Preprocessing Steps Performed

# See function "preproc" in funs.R and dataPrep.R
  1. Verbal fluency scores normalized for each version (K,P,R,S) by subtracting mean per version from each score.
  2. Concatenation of Wellbeing and Training Motivation and removal of entries where both are 50, assuming that they might have been chosen out of “laziness”.
  3. Concatenation of Operation Span & N-Back accuracy scores due to high correlation.
  4. Logarithmising all RTs.
  5. Z-normalization of all cognitive variables, so that ranges are comparable.
  6. For the overall models, RTs were multiplied by -1, so that the scaling is comparable between accuracy and RTs (higher = better).
  7. The very first time point was removed as familiarization trial (in which participants also received feedback from the investigators).




Overview of Variables After Preprocessing


Data Description

(dataframe df)

  • userCode: subject ID
  • age: Subject’s age in years at study commencement.
  • educ_y: Number of years in education up to university. Advanced trainings included.
  • Act_total: Number of hobbies practiced on a regular basis as assessed via questionnaire.
  • BLP: Score for multilingualism.
  • Time: number of measurement per subject
  • Group: Experimental group. PASV = Passive control group (social exchange only); ACTV= Active control group (strategy games); LANG= Experimental group (Spanish training).
  • isLANG: Binary factor of whether or not participants were part of the language group or not (1 = yes)
  • isLANG0: Ordered factor of whether or not participants were part of the language group or not (1 = yes)
  • isACTV: Binary factor of whether or not participants were part of the gaming group or not (1 = yes)
  • isACTV0: Ordered factor of whether or not participants were part of the gaming group or not (1 = yes)
  • isPASV: Binary factor of whether or not participants were part of the movies group or not (1 = yes)
  • isPASV0: Ordered factor of whether or not participants were part of the movies group or not (1 = yes)
  • class: Groups of 5-12 individuals, split up by experimental condition (LANG: 4 groups, PASV: 2 groups, ACTV: 2 groups).
  • gend: Gender as binary (1 = female)
  • socioAffect: Mean score 1-100 of overall wellbeing and training motivation during the week in question (100 is highest).
  • sRWT: Number of words produced within 1 minute. Score for verbal fluency (z-score).
  • sWM: Sensitivity/accuracy score over both working memory tasks (Operation Span & 2-Back)(z-score).
  • rtWM: Logarithmised reaction time for 2-Back task (z-score).
  • sDivAtt: Sensitivity/accuracy score in Divided Attention task. Not included in analysis due to pronounced ceiling effects (z-score).
  • rtDivAtt Logarithmised reaction time for Divided Attention task (z-score).
  • rtAlert: Logarithmised mean reaction time for Alertness task (z-score).
  • start.event: Variable necessary for autocorrelation in GAMMs.
  • initLevel_sWM: Baseline Working Memory Accuracy from GAMs per subject
  • initLevel_rtWM: Baseline Working Memory RT from GAMs per subject
  • initLevel_rtDivAtt: Baseline Divided Attention RT from GAMs per subject
  • initLevel_rtAlert: Baseline Alertness RT from GAMs per subject
  • initLevel_sRWT: Baseline Verbal Fluency from GAMs per subject






Data Excerpt

##   userCode age educ_y Act_total     BLP Time Group isLANG isLANG0 isACTV
## 1    0uh8j  67     12        14 20.9975    0  PASV      0       0      0
## 2    0uh8j  67     12        14 20.9975    1  PASV      0       0      0
## 3    0uh8j  67     12        14 20.9975    2  PASV      0       0      0
## 4    0uh8j  67     12        14 20.9975    3  PASV      0       0      0
## 5    0uh8j  67     12        14 20.9975    4  PASV      0       0      0
## 6    0uh8j  67     12        14 20.9975    5  PASV      0       0      0
##   isACTV0 isPASV isPASV0 class gend socioAffect     sRWT       sWM        rtWM
## 1       0      1       1  Fri4    0      88.150 9.192869 0.7929436 -0.18233717
## 2       0      1       1  Fri4    0      84.720 8.098263 0.6200322 -0.22773324
## 3       0      1       1  Fri4    0      79.360 9.165966 0.7113189 -0.04191973
## 4       0      1       1  Fri4    0      77.345 9.098263 0.6885683 -0.06018706
## 5       0      1       1  Fri4    0      63.410 9.510919 0.7384281  0.01091307
## 6       0      1       1  Fri4    0      76.890 8.192869 0.7502542 -0.19012933
##    sDivAtt  rtDivAtt   rtAlert start.event initLevel_sWM initLevel_rtWM
## 1 3.213587 0.4569827 -1.199853        TRUE     0.7240647     -0.1335443
## 2 1.969167 0.3838274 -1.068783       FALSE     0.7240647     -0.1335443
## 3 4.762174 0.4402467 -1.271873       FALSE     0.7240647     -0.1335443
## 4 2.036882 0.5306744 -1.299131       FALSE     0.7240647     -0.1335443
## 5 4.074072 0.3911455 -1.310328       FALSE     0.7240647     -0.1335443
## 6 3.494582 0.2529916 -1.294960       FALSE     0.7240647     -0.1335443
##   initLevel_rtDivAtt initLevel_rtAlert initLevel_sRWT
## 1          0.4469415          -1.16499       8.798443
## 2          0.4469415          -1.16499       8.798443
## 3          0.4469415          -1.16499       8.798443
## 4          0.4469415          -1.16499       8.798443
## 5          0.4469415          -1.16499       8.798443
## 6          0.4469415          -1.16499       8.798443

Descriptive Statistics

##    userCode              age            educ_y        Act_total   
##  Length:1753        Min.   :64.00   Min.   : 8.00   Min.   :10.0  
##  Class :character   1st Qu.:66.00   1st Qu.:10.00   1st Qu.:15.0  
##  Mode  :character   Median :67.00   Median :12.50   Median :17.0  
##                     Mean   :68.42   Mean   :12.65   Mean   :16.4  
##                     3rd Qu.:70.00   3rd Qu.:15.00   3rd Qu.:19.0  
##                     Max.   :75.00   Max.   :20.00   Max.   :22.0  
##                                                                   
##       BLP               Time        Group         isLANG       isLANG0
##  Min.   : 0.1135   Min.   : 0.00   LANG:809   Min.   :0.0000   0:944  
##  1st Qu.:14.7550   1st Qu.: 7.00   ACTV:489   1st Qu.:0.0000   1:809  
##  Median :26.3095   Median :14.00   PASV:455   Median :0.0000          
##  Mean   :28.8126   Mean   :13.88              Mean   :0.4615          
##  3rd Qu.:40.2480   3rd Qu.:21.00              3rd Qu.:1.0000          
##  Max.   :90.0965   Max.   :28.00              Max.   :1.0000          
##                                                                       
##      isACTV      isACTV0      isPASV       isPASV0     class          
##  Min.   :0.000   0:1264   Min.   :0.0000   0:1298   Length:1753       
##  1st Qu.:0.000   1: 489   1st Qu.:0.0000   1: 455   Class :character  
##  Median :0.000            Median :0.0000            Mode  :character  
##  Mean   :0.279            Mean   :0.2596                              
##  3rd Qu.:1.000            3rd Qu.:1.0000                              
##  Max.   :1.000            Max.   :1.0000                              
##                                                                       
##       gend         socioAffect          sRWT             sWM        
##  Min.   :0.0000   Min.   :  8.72   Min.   : 4.098   Min.   :0.1573  
##  1st Qu.:0.0000   1st Qu.: 63.17   1st Qu.:10.193   1st Qu.:0.6856  
##  Median :0.0000   Median : 79.17   Median :12.511   Median :0.7806  
##  Mean   :0.4261   Mean   : 75.82   Mean   :12.663   Mean   :0.7691  
##  3rd Qu.:1.0000   3rd Qu.: 91.26   3rd Qu.:15.098   3rd Qu.:0.8703  
##  Max.   :1.0000   Max.   :100.00   Max.   :22.511   Max.   :1.0000  
##                   NA's   :71       NA's   :8        NA's   :4       
##       rtWM            sDivAtt          rtDivAtt            rtAlert       
##  Min.   :-1.0221   Min.   :0.2281   Min.   :-0.622314   Min.   :-1.4578  
##  1st Qu.:-0.5757   1st Qu.:5.8081   1st Qu.:-0.225384   1st Qu.:-1.2609  
##  Median :-0.4072   Median :6.9844   Median : 0.005215   Median :-1.1780  
##  Mean   :-0.3848   Mean   :6.7427   Mean   : 0.000403   Mean   :-1.1649  
##  3rd Qu.:-0.2199   3rd Qu.:8.0824   3rd Qu.: 0.231475   3rd Qu.:-1.0916  
##  Max.   : 0.4894   Max.   :8.2217   Max.   : 0.555148   Max.   :-0.2886  
##  NA's   :12        NA's   :16       NA's   :7           NA's   :2        
##  start.event     initLevel_sWM    initLevel_rtWM    initLevel_rtDivAtt
##  Mode :logical   Min.   :0.2320   Min.   :-0.7742   Min.   :-0.27715  
##  FALSE:1692      1st Qu.:0.5970   1st Qu.:-0.5167   1st Qu.:-0.02521  
##  TRUE :61        Median :0.7045   Median :-0.3222   Median : 0.12078  
##                  Mean   :0.6735   Mean   :-0.3086   Mean   : 0.11026  
##                  3rd Qu.:0.7717   3rd Qu.:-0.1399   3rd Qu.: 0.23448  
##                  Max.   :0.9344   Max.   : 0.3750   Max.   : 0.44694  
##                                                                       
##  initLevel_rtAlert initLevel_sRWT  
##  Min.   :-1.3962   Min.   : 6.217  
##  1st Qu.:-1.2154   1st Qu.: 9.326  
##  Median :-1.1337   Median :11.498  
##  Mean   :-1.1332   Mean   :11.391  
##  3rd Qu.:-1.0510   3rd Qu.:13.089  
##  Max.   :-0.6541   Max.   :17.890  
##

RQ1 & RQ2: Group Differences in Cognitive Development

1. What is the impact of training type (LANG, ACTV or PASV) on dynamic processes of change in cognitive functioning? 2. Is there a difference between the three training types (LANG, ACTV or PASV) in terms of the overall trend of change through time?



M2: Overall Training Effect per Group

The GAMM model to test the group effect on overall cognitive performance was computed as follows:

m2_old = bam(score ~ 
            isACTV0 +
            isPASV0 +
            s(Time) +
            s(Time, by = isACTV0) +
            s(Time, by = isPASV0) +
            s(Time, userCode, bs = "fs", m = 1) + 
            s(Time, Tasks, bs = "fs", m = 1),
          data = df_overall, family = "scat", discrete = T, nthreads = 2)

r1 = start_value_rho(m2_old)

m2 = bam(score ~ 
            isACTV0 +
            isPASV0 +
            s(Time) +
            s(Time, by = isACTV0) +
            s(Time, by = isPASV0) +
            s(Time, userCode, bs = "fs", m = 1) + 
            s(Time, Tasks, bs = "fs", m = 1),
          data = df_overall, rho = r1,  AR.start = df_overall$start.event,
          family = "scat", discrete = T, nthreads = 2)

# The code is executed in `3_GAMMs.R`.

Explanation:

  • df_overallincludes all tasks per time and subject in long format.
  • isACTV0 and isPASV0 are parametric terms that describe how the group’s intercepts differs from that of the reference group (LANG).
  • s(Time) = The overall training effect independent of groups (i.e. the average practice effect).
  • s(Time, by = isACTV0) and s(Time, by = isPASV0) = The difference smooths for the effect of time, i.e. how the training effect in ACTV and PASV was different from the overall training gain in LANG.
  • s(Time, userCode, bs = "fs", m = 1) = Factor smooth interaction (i.e. random effects) for the effect of time within a given subject at the reference level (LANG).
  • s(Time, Tasks, bs = "fs", m = 1) = Random smooths for task type (i.e. a separate smooth for each of the 5 tasks).
  • rho and AR.start remove autocorrelation based on the rho value of the first model.
  • family = "scat" applies the scaled-t model family, which yielded a more normal distribution of residuals.



M2: Visualization of Overall Trajectories



See 3_GAMMs.R for visualization code.



M2: Output of Overall Model


Family: gaussian 
Link function: identity 

Formula:
score ~ isACTV0 + isPASV0 + s(Time) + s(Time, by = isACTV0) + 
    s(Time, by = isPASV0) + s(Time, userCode, bs = "fs", m = 1)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.11757    0.09728   1.209   0.2269  
isACTV01    -0.16072    0.15816  -1.016   0.3096  
isPASV01    -0.35316    0.16131  -2.189   0.0286 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                    edf  Ref.df      F p-value    
s(Time)           3.594   4.585 17.192  <2e-16 ***
s(Time):isACTV01  1.000   1.000  0.014   0.907    
s(Time):isPASV01  1.357   1.632  0.209   0.679    
s(Time,userCode) 53.967 547.000  1.413  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.302   Deviance explained = 30.6%
fREML = 8357.9  Scale est. = 0.64752   n = 8732

There is no significant difference in the overall cognitive trajectories between LANG and PASV or between LANG and ACTV.



M2: Model Criticism for Overall Model

load(paste(outputFolder, "tables_and_summaries", "m2.rda", sep = "/"))

par(mfrow = c(1,2))
  acf_resid(m2, main = "ACT Plot of Model Residuals")
  qqp(resid(m2), id = F, main = "Quantile-Comparison Plot")


There is no noticeable autocorrelation and the model residuals, and the residuals are normally distributed.



M2: Predicting Group Differences at T=0

langPred = get_predictions(m2, cond = list(isACTV0 = 0, 
                                       isPASV0 = 0,
                                       Time = 0),
                           rm.ranef = T, print.summary = F)

actvPred = get_predictions(m2, cond = list(isACTV0 = 1, 
                                       isPASV0 = 0,
                                       Time = 0),
                           rm.ranef = T, print.summary = F)

pasvPred = get_predictions(m2, cond = list(isACTV0 = 0, 
                                       isPASV0 = 1,
                                       Time = 0),
                           rm.ranef = T, print.summary = F)

At time point zero, overall cognitive performance at the reference level was -0.3, for group ACTV it was -0.47 and for PASV it was -0.71.

M3: Training Effect per Cognitive Task

The GAMM model to test the group effect on each cognitive measure was computed as follows:

# Model that is not corrected for autocorrelation
m3_old = bam(score ~ 
               
               # smooths for reference level LANG
               s(Time, by = Tasks) + 
               Tasks + 
               
               # smooths for ACTV as ordered factors
               s(Time, by = isACTVissWM0) +
               isACTVissWM0 + 
               s(Time, by = isACTVisrtWM0) +
               isACTVisrtWM0 + 
               s(Time, by = isACTVisrtDivAtt0) +
               isACTVisrtDivAtt0 + 
               s(Time, by = isACTVisrtAlert0) +
               isACTVisrtAlert0 + 
               s(Time, by = isACTVissRWT0) +
               isACTVissRWT0 + 
               
               # smooths for PASV as ordered factors
               s(Time, by = isPASVissWM0) +
               isPASVissWM0 + 
               s(Time, by = isPASVisrtWM0) +
               isPASVisrtWM0 + 
               s(Time, by = isPASVisrtDivAtt0) +
               isPASVisrtDivAtt0 + 
               s(Time, by = isPASVisrtAlert0) +
               isPASVisrtAlert0 + 
               s(Time, by = isPASVissRWT0) +
               isPASVissRWT0 +
               
               # random effects for subject and task in LANG
               s(Time, userCode, by = Tasks, bs = "fs", m = 1),
             data = df_overall, family = "scat", discrete = T, nthreads = 2)

r1 = start_value_rho(m3_old)

# Model corrected for autocorrelation
m3 = bam(score ~ 
           # smooths for reference level LANG
           s(Time, by = Tasks) + 
           Tasks + 
           
           # smooths for ACTV as ordered factors
           s(Time, by = isACTVissWM0) +
           isACTVissWM0 + 
           s(Time, by = isACTVisrtWM0) +
           isACTVisrtWM0 + 
           s(Time, by = isACTVisrtDivAtt0) +
           isACTVisrtDivAtt0 + 
           s(Time, by = isACTVisrtAlert0) +
           isACTVisrtAlert0 + 
           s(Time, by = isACTVissRWT0) +
           isACTVissRWT0 + 
           
           # smooths for PASV as ordered factors
           s(Time, by = isPASVissWM0) +
           isPASVissWM0 + 
           s(Time, by = isPASVisrtWM0) +
           isPASVisrtWM0 + 
           s(Time, by = isPASVisrtDivAtt0) +
           isPASVisrtDivAtt0 + 
           s(Time, by = isPASVisrtAlert0) +
           isPASVisrtAlert0 + 
           s(Time, by = isPASVissRWT0) +
           isPASVissRWT0 +
           
           # random effects for subject and task in LANG
           s(Time, userCode, by = Tasks, bs = "fs", m = 1),
         data = df_overall, family = "scat", rho = r1,  AR.start = df_overall$start.event,
         discrete = T, nthreads = 2)


# The code is executed in GAMM.R.

Explanation:

  • The term Tasks models a separate intercept for each of the cognitive measures.
  • s(Time, by = Tasks) models a smooth for the overall training/practice effect per cognitive measure, independent of the training type.
  • s(Time, by = isACTVissWM0) and similar terms model a difference smooth per group (using ordered factors) and task type, while isACTVissWM0 models the respective intercept difference per group and task. For example, this smooth models how the trajectory of Working Memory Accuracy over time differed between group ACTV and group LANG.
  • All other term explanations are the same as for the overall model.



M3: Output from Model per Cognitive Task


Family: gaussian 
Link function: identity 

Formula:
score ~ s(Time, by = Tasks) + Tasks + s(Time, by = isACTVissWM0) + 
    isACTVissWM0 + s(Time, by = isACTVisrtWM0) + isACTVisrtWM0 + 
    s(Time, by = isACTVisrtDivAtt0) + isACTVisrtDivAtt0 + s(Time, 
    by = isACTVisrtAlert0) + isACTVisrtAlert0 + s(Time, by = isACTVissRWT0) + 
    isACTVissRWT0 + s(Time, by = isPASVissWM0) + isPASVissWM0 + 
    s(Time, by = isPASVisrtWM0) + isPASVisrtWM0 + s(Time, by = isPASVisrtDivAtt0) + 
    isPASVisrtDivAtt0 + s(Time, by = isPASVisrtAlert0) + isPASVisrtAlert0 + 
    s(Time, by = isPASVissRWT0) + isPASVissRWT0 + s(Time, userCode, 
    bs = "fs", m = 1)

Parametric coefficients:
                    Estimate Std. Error t value Pr(>|t|)   
(Intercept)         0.179656   0.109587   1.639  0.10117   
TasksrtWM           0.004251   0.079845   0.053  0.95754   
TasksrtDivAtt      -0.126477   0.079953  -1.582  0.11371   
TasksrtAlert       -0.114533   0.079623  -1.438  0.15035   
TaskssRWT          -0.090288   0.079843  -1.131  0.25816   
isACTVissWM01      -0.225489   0.178000  -1.267  0.20526   
isACTVisrtWM01     -0.430346   0.178033  -2.417  0.01566 * 
isACTVisrtDivAtt01  0.065242   0.178001   0.367  0.71398   
isACTVisrtAlert01  -0.081396   0.178115  -0.457  0.64769   
isACTVissRWT01     -0.119142   0.177982  -0.669  0.50326   
isPASVissWM01      -0.544538   0.181571  -2.999  0.00272 **
isPASVisrtWM01     -0.291410   0.181616  -1.605  0.10863   
isPASVisrtDivAtt01 -0.393329   0.181495  -2.167  0.03025 * 
isPASVisrtAlert01  -0.226900   0.181731  -1.249  0.21186   
isPASVissRWT01     -0.288972   0.181804  -1.589  0.11199   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                              edf  Ref.df      F  p-value    
s(Time):TaskssWM            3.182   4.062 12.076  < 2e-16 ***
s(Time):TasksrtWM           1.001   1.001  6.446 0.011100 *  
s(Time):TasksrtDivAtt       2.418   3.078  6.751 0.000137 ***
s(Time):TasksrtAlert        1.001   1.002  0.600 0.438244    
s(Time):TaskssRWT           2.029   2.570  6.216 0.001130 ** 
s(Time):isACTVissWM01       1.000   1.000  0.122 0.726602    
s(Time):isACTVisrtWM01      1.000   1.000  0.014 0.905549    
s(Time):isACTVisrtDivAtt01  1.000   1.000  0.018 0.892740    
s(Time):isACTVisrtAlert01   1.821   2.302  0.763 0.469820    
s(Time):isACTVissRWT01      1.000   1.000  0.082 0.774090    
s(Time):isPASVissWM01       1.000   1.000  1.154 0.282768    
s(Time):isPASVisrtWM01      1.000   1.000  0.001 0.973245    
s(Time):isPASVisrtDivAtt01  1.000   1.000  0.137 0.711491    
s(Time):isPASVisrtAlert01   2.107   2.679  1.151 0.340404    
s(Time):isPASVissRWT01      1.685   2.109  0.560 0.571820    
s(Time,userCode)           54.105 546.000  1.469  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.315   Deviance explained = 32.2%
fREML = 8375.6  Scale est. = 0.63688   n = 8732
# Code for summary above was executed in 3_GAMMs.R

M3: Visualization per Task



The visualization code can be found in 3_GAMMs.R



M3: Difference Plots per Task

(Note: Values are standardized, such that higher scores represent better performance, also for RTs.)

par(mfrow = c(2,5))
    plot_diff(m3, 
              view = "Time", 
              comp = list(isPASVissWM0 = c("0", "1")),
              cond = list(Tasks = "sWM"),
              ylab = "LANG minus PASV",
              main = "Working Memory Acc.",
              rm.ranef = T,
              print.summary = F)
  
      plot_diff(m3, 
              view = "Time", 
              comp = list(isPASVisrtWM0 = c("0", "1")),
              cond = list(Tasks = "rtWM"),
              ylab = "",
              main = "Working Memory RT",
              rm.ranef = T,
              print.summary = F)
  
      plot_diff(m3, 
              view = "Time", 
              comp = list(isPASVisrtDivAtt0 = c("0", "1")),
              cond = list(Tasks = "rtDivAtt"),
              ylab = "",
              main = "Divided Attention RT",
              rm.ranef = T,
              print.summary = F)
      
      plot_diff(m3, 
              view = "Time", 
              comp = list(isPASVisrtAlert0 = c("0", "1")),
              cond = list(Tasks = "rtAlert"),
              ylab = "",
              main = "Alertness RT",
              rm.ranef = T,
              print.summary = F)
      
      plot_diff(m3, 
              view = "Time", 
              comp = list(isPASVissRWT0 = c("0", "1")),
              cond = list(Tasks = "sRWT"),
              ylab = "",
              main = "Verbal Fluency",
              rm.ranef = T,
              print.summary = F)
      
      plot_diff(m3, 
              view = "Time", 
              comp = list(isACTVissWM0 = c("0", "1")),
              cond = list(Tasks = "sWM"),
              ylab = "LANG minus ACTV",
              main = "",
              rm.ranef = T,
              print.summary = F)
      
      plot_diff(m3, 
              view = "Time", 
              comp = list(isACTVisrtWM0 = c("0", "1")),
              cond = list(Tasks = "rtWM"),
              ylab = "",
              main = "",
              rm.ranef = T,
              print.summary = F)
      
      plot_diff(m3, 
              view = "Time", 
              comp = list(isACTVisrtDivAtt0 = c("0", "1")),
              cond = list(Tasks = "rtDivAtt"),
              ylab = "",
              main = "",
              rm.ranef = T,
              print.summary = F)
      
      plot_diff(m3, 
              view = "Time", 
              comp = list(isACTVisrtAlert0 = c("0", "1")),
              cond = list(Tasks = "rtAlert"),
              ylab = "",
              main = "",
              rm.ranef = T,
              print.summary = F)
      
    plot_diff(m3, 
              view = "Time", 
              comp = list(isACTVissRWT0 = c("0", "1")),
              cond = list(Tasks = "sRWT"),
              ylab = "",
              main = "",
              rm.ranef = T,
              print.summary = F)


Result:

  • Working Memory Accuracy was higher in group LANG than in group PASV throughout the training (= selection bias), but the difference kept decreasing.
  • Divided Attention RT was significantly better (i.e. lower) in group LANG compared to PASV after week 7, but here, too, there was a selection bias as shown in differences at T=0.
  • Group LANG also outperformed group ACTV in terms of Working Memory RT and significantly so after week 3, but again, values at T=0 were already higher at the beginning, pointing to a selection bias.



M3: Model Criticism

par(mfrow = c(1,2))
  acf_resid(m3, main = "ACT Plot of Model Residuals")
  qqp(resid(m3), id = F, main = "Quantile-Comparison Plot")


No autocorrelation. Residuals normally distributed.




M2 & M3 Results Summarized

  • Overall (M3)
    • Significant practice effect.
    • Group LANG had higher intercept than PASV.
    • No overall group difference in smooth over time (i.e. no difference in cognitive development between ACTV and LANG or between PASV and LANG).
  • Working Memory Accuracy (M2)
    • Significant overall practice effect.
    • Group PASV had significantly lower intercept than LANG.
    • No significant difference in smooth over time between LANG and PASV or between LANG and ACTV.
    • Difference plot shows significantly higher accuracy in group LANG compared to PASV throughout the training when accounting for uncertainty in intercept.
  • Working Memory RT (M2)
    • Significant overall practice effect.
    • Group LANG had a significantly higher intercept than group ACTV.
    • No group difference in smooth over time.
  • Divided Attention RT (M2)
    • Significant overall practice effect.
    • Group LANG had significantly higher intercept than group PASV.
    • No group difference in smooth over time.
    • Difference plot shows significantly better performance in group LANG compared to PASV after week 7 when accounting for uncertainty in intercept.
  • Alertness Reaction Time (M2)
    • No significant practice effect.
    • No group difference in intercept.
    • No group difference in smooth over time.
  • Verbal Fluency (M2)
    • Significant practice effect.
    • No group difference in intercept.
    • No group difference in smooth over time.










RQ3 & RQ4: Individual Differences in Cognitive Benefit

Ideally, groups should have been identical in terms of ID variables, cognitive baseline and socio-affect in order to be able to objectively quantify the transfer effect of LANG and ACTV training modalities. While they were comparable in ID variables (i.e. no selection bias) and socio-affect, there were slight differences in cognitive baseline (LANG starting with higher baseline). Thus, in order to assess the pure effect of the LANG / ACTV trainings, a model was computed in which baseline level, time and their interaction were inserted as predictors.



M4_te: Overall Interaction Time x Baseline Level Using te()

The following was the model code for M4_te:

m4_te_old = bam(score ~ 
                  te(Time, initLevel) +
                  te(Time, initLevel, by = isACTV) +
                  te(Time, initLevel, by = isPASV) +
                s(Time, userCode, bs = "fs", m = 1),
              data = df_overall, family = "scat", discrete = T, nthreads = 7)

r1 = start_value_rho(m4_te_old)

m4_te = bam(score ~ 
                  te(Time, initLevel) +
                  te(Time, initLevel, by = isACTV) +
                  te(Time, initLevel, by = isPASV) +
              s(Time, userCode, bs = "fs", m = 1),
          data = df_overall, rho = r1,  AR.start = df_overall$start.event,
          family = "scat", discrete = T, nthreads = 7)

# Model run in 3_GAMMs.R


Explanation:

te() produces a full tensor product interaction, i.e. an interaction that includes the main effects of the respective terms.



M4_te: Model Output


Family: gaussian 
Link function: identity 

Formula:
score ~ te(Time, initLevel) + te(Time, initLevel, by = isACTV) + 
    te(Time, initLevel, by = isPASV) + s(Time, userCode, bs = "fs", 
    m = 1)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.04110    0.04833    0.85    0.395

Approximate significance of smooth terms:
                             edf  Ref.df      F  p-value    
te(Time,initLevel)        15.994  18.393 79.736  < 2e-16 ***
te(Time,initLevel):isACTV  9.393   9.850  3.607 5.24e-05 ***
te(Time,initLevel):isPASV  6.491   7.429  2.724  0.00674 ** 
s(Time,userCode)          86.508 546.000  1.037  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.595   Deviance explained =   60%
fREML = 7659.1  Scale est. = 0.39445   n = 8732



par(mfrow = c(1,2))
  
  # te(Time,initLevel):isPASV
  pvisgam(m4_te, select = 3,
       main = "Smooth over Baseline Level \n in PASV",
       view = c("Time", "initLevel"),
       zlim = c(-1.15, 1.15),
       print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).
  # te(Time,initLevel):isACTV
  pvisgam(m4_te, select = 2,
       main = "Smooth over Baseline Level \n in ACTV",
       view = c("Time", "initLevel"),
       zlim = c(-2.23, 2.23),
       print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).

The effect of baseline level and its interaction with time was significantly different in groups ACTV and PASV from LANG such that individuals with low cognitive baseline performance improved more in group LANG than in both ACTV and PASV. However, it is possible that the model was largely influenced by outliers on the lower spectrum of cognitive baselines, especially the two outliers where baseline was < 3.5SD.

See the distribution of baseline levels:

dat = df_overall %>%
  group_by(userCode, Tasks) %>%
  summarise(initLevel = mean(initLevel))
hist(dat$initLevel, 50, xlim = c(-4,4),xlab = "Baseline Level", main = "Histogram of Baseline Levels")



M4_threeOut (without outliers < -3.5SD)

Model M4 was run again without outliers < -3.5SD.


Family: gaussian 
Link function: identity 

Formula:
score ~ te(Time, initLevel) + te(Time, initLevel, by = isACTV) + 
    te(Time, initLevel, by = isPASV) + s(Time, userCode, bs = "fs", 
    m = 1)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.04309    0.04830   0.892    0.372

Approximate significance of smooth terms:
                             edf  Ref.df      F  p-value    
te(Time,initLevel)        16.179  18.690 78.355  < 2e-16 ***
te(Time,initLevel):isACTV  9.424   9.864  3.787 2.65e-05 ***
te(Time,initLevel):isPASV  6.917   7.912  1.996   0.0399 *  
s(Time,userCode)          85.599 547.000  1.029  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.587   Deviance explained = 59.2%
fREML = 7561.5  Scale est. = 0.39117   n = 8678
par(mfrow = c(1,2))
  
  # te(Time,initLevel):isPASV
  pvisgam(m4_threeOut, select = 3,
       main = "Smooth over Baseline Level \n in PASV",
       view = c("Time", "initLevel"),
       zlim = c(-0.95, 0.95),
       print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).
  # te(Time,initLevel):isACTV
  pvisgam(m4_threeOut, select = 2,
       main = "Smooth over Baseline Level \n in ACTV",
       view = c("Time", "initLevel"),
       zlim = c(-1.00, 1.00),
       print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).



Model Comparison M2 (Time only) vs. M4 (with Baseline Cognition)

M2 and M4 were modeled again assuming all effects to be random effects by using select = T (see 3_GAMMs.R).

The models used binary curves for group differences and excluded outliers < -3.5SD.

Model comparison (using compareML()) results in the following values:

load(file = paste(outputFolder, "tables_and_summaries", "m2_threeout_alt.rda", sep = "/"))
load(file = paste(outputFolder, "tables_and_summaries", "m4_threeout_alt.rda", sep = "/"))

comparison = compareML(m2.alt, m4.alt)
## m2.alt: score ~ isACTV0 + isPASV0 + s(Time) + s(Time, by = isACTV0) + 
##     s(Time, by = isPASV0) + s(Time, userCode, bs = "fs", m = 1)
## 
## m4.alt: score ~ te(Time, initLevel) + te(Time, initLevel, by = isACTV) + 
##     te(Time, initLevel, by = isPASV) + s(Time, userCode, bs = "fs", 
##     m = 1)
## 
## Chi-square test of fREML scores
## -----
##    Model    Score Edf Difference    Df  p.value Sig.
## 1 m2.alt 8235.667  11                               
## 2 m4.alt 7575.740  12    659.927 1.000  < 2e-16  ***
## 
## AIC difference: 1352.91, model m4.alt has lower AIC.

M4 had a lower AIC and fREML score, and is therefore preferred.




M5_threeOut: Adding Individual Tasks to M4_threeOut

To investigate the effect of baseline level per cognitive task, the following was the model code for M5_threeOut. The outliers from model M4_threeOut were removed here, too.

m5_threeOut_old = bam(score ~ 
                      te(Time, initLevel, by = Tasks) +
                      te(Time, initLevel, by = isACTVissWM) +
                      te(Time, initLevel, by = isACTVisrtWM) +
                      te(Time, initLevel, by = isACTVisrtDivAtt) +
                      te(Time, initLevel, by = isACTVisrtAlert) +
                      te(Time, initLevel, by = isACTVissRWT) +
                      
                      te(Time, initLevel, by = isPASVissWM) +
                      te(Time, initLevel, by = isPASVisrtWM) +
                      te(Time, initLevel, by = isPASVisrtDivAtt) +
                      te(Time, initLevel, by = isPASVisrtAlert) +
                      te(Time, initLevel, by = isPASVissRWT) +
                      
                      s(Time, userCode, bs = "fs", m = 1),
                    data = df_overall3, family = "scat", discrete = T, nthreads = 7)

r1 = start_value_rho(m5_threeOut_old)

m5_threeOut = bam(score ~ 
                  te(Time, initLevel, by = Tasks) +
                  te(Time, initLevel, by = isACTVissWM) +
                  te(Time, initLevel, by = isACTVisrtWM) +
                  te(Time, initLevel, by = isACTVisrtDivAtt) +
                  te(Time, initLevel, by = isACTVisrtAlert) +
                  te(Time, initLevel, by = isACTVissRWT) +
                  
                  te(Time, initLevel, by = isPASVissWM) +
                  te(Time, initLevel, by = isPASVisrtWM) +
                  te(Time, initLevel, by = isPASVisrtDivAtt) +
                  te(Time, initLevel, by = isPASVisrtAlert) +
                  te(Time, initLevel, by = isPASVissRWT) +
                  s(Time, userCode, bs = "fs", m = 1),
                data = df_overall3, rho = r1,  AR.start = df_overall3$start.event,
                family = "scat", discrete = T, nthreads = 7)




M5_threeOut: Model Output


Family: gaussian 
Link function: identity 

Formula:
score ~ te(Time, initLevel, by = Tasks) + te(Time, initLevel, 
    by = isACTVissWM) + te(Time, initLevel, by = isACTVisrtWM) + 
    te(Time, initLevel, by = isACTVisrtDivAtt) + te(Time, initLevel, 
    by = isACTVisrtAlert) + te(Time, initLevel, by = isACTVissRWT) + 
    te(Time, initLevel, by = isPASVissWM) + te(Time, initLevel, 
    by = isPASVisrtWM) + te(Time, initLevel, by = isPASVisrtDivAtt) + 
    te(Time, initLevel, by = isPASVisrtAlert) + te(Time, initLevel, 
    by = isPASVissRWT) + s(Time, userCode, bs = "fs", m = 1)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.02129    0.05147   0.414    0.679

Approximate significance of smooth terms:
                                       edf  Ref.df      F  p-value    
te(Time,initLevel):TaskssWM         14.032  16.999 19.368  < 2e-16 ***
te(Time,initLevel):TasksrtWM         5.410   6.359 61.125  < 2e-16 ***
te(Time,initLevel):TasksrtDivAtt    10.851  13.293 45.408  < 2e-16 ***
te(Time,initLevel):TasksrtAlert     11.639  13.477 38.071  < 2e-16 ***
te(Time,initLevel):TaskssRWT         7.779  10.118 39.590  < 2e-16 ***
te(Time,initLevel):isACTVissWM       4.000   4.001  2.924 0.019795 *  
te(Time,initLevel):isACTVisrtWM      9.003   9.633  4.638 4.39e-07 ***
te(Time,initLevel):isACTVisrtDivAtt  4.000   4.000  0.380 0.823103    
te(Time,initLevel):isACTVisrtAlert   6.990   7.845  2.468 0.012255 *  
te(Time,initLevel):isACTVissRWT      7.996   8.922  1.611 0.113826    
te(Time,initLevel):isPASVissWM       6.009   6.702  3.413 0.001668 ** 
te(Time,initLevel):isPASVisrtWM      9.324   9.832  4.463 1.00e-06 ***
te(Time,initLevel):isPASVisrtDivAtt  8.117   9.026  1.985 0.033403 *  
te(Time,initLevel):isPASVisrtAlert   7.959   9.186  2.253 0.014037 *  
te(Time,initLevel):isPASVissRWT      8.009   8.879  3.341 0.000477 ***
s(Time,userCode)                    96.699 546.000  1.329  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =   0.63   Deviance explained = 63.9%
fREML = 7431.8  Scale est. = 0.35075   n = 8678
par(mfrow = c(2,5))

  # te(Time,initLevel):isPASVissWM 
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 11,
          main = "Working Memory Acc. \n in PASV",
          ylab = "Baseline Level",
          zlim = c(-1.7, 1.7),
          print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).
  # te(Time,initLevel):isPASVisrtWM 
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 12,
          main = "Working Memory RT \n in PASV",
          ylab = "Baseline Level",
          zlim = c(-7.62, 7.62),
          print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).
  # te(Time,initLevel):isPASVisrtDivAtt
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 13,
          main = "Divided Attention RT \n in PASV",
          ylab = "Baseline Level",
          zlim = c(-2.54, 2.54),
          print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).
  # te(Time,initLevel):isPASVisrtAlert 
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 14,
          main = "Alertness RT \n in PASV",
          ylab = "Baseline Level",
          zlim = c(-2.18, 2.18),
          print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).
    # te(Time,initLevel):isPASVissRWT
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 15,
          main = "Verbal Fluency \n in PASV",
          ylab = "Baseline Level",
          zlim = c(-1.98, 1.98),
          print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).
  # te(Time,initLevel):isACTVissWM 
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 6,
          main = "Working Memory Acc. \n in ACTV",
          ylab = "Baseline Level",
          zlim = c(-0.70, 0.70),
          print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).
  # te(Time,initLevel):isACTVisrtWM 
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 7,
          main = "Working Memory RT \n in ACTV",
          ylab = "Baseline Level",
          zlim = c(-2.6, 2.6),
          print.summary = F)
  
  # te(Time,initLevel):isACTVisrtDivAtt
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 8,
          main = "Divided Attention RT \n in ACTV",
          ylab = "Baseline Level",
          zlim = c(-2.54, 2.54),
          print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).
  # te(Time,initLevel):isACTVisrtAlert 
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 9,
          main = "Alertness RT \n in ACTV",
          ylab = "Baseline Level",
          zlim = c(-2.00, 2.00),
          print.summary = F)
  
  # te(Time,initLevel):isACTVissRWT 
  pvisgam(m5_threeOut, view = c("Time", "initLevel"),
          select = 10,
          main = "Verbal Fluency \n in ACTV",
          ylab = "Baseline Level",
          zlim = c(-2.9, 2.9),
          print.summary = F)
## Warning in gradientLegend(round(c(min.z, max.z), 3), n.seg = 3, pos = 0.875, :
## Increase right margin to fit labels or decrease the number of decimals, see
## help(gradientLegend).